The response to dynamical modulation of the optical lattice for fermions in the 

Hubbard model 
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Fermionic atoms in a periodic optical lattice provide a realization of the single-band Hubbard 
model. Using Quantum Monte Carlo simulations along with the Maximum Entropy Method, we 
evaluate the effect of a time-dependent perturbative modulation of the optical lattice amplitude 
on atomic correlations, revealed in the fraction of doubly-occupied sites. Our treatment extends 
previous approaches which neglected the time dependence of the on-site interaction, and shows 
that this term changes the results in a quantitatively significant way. The effect of modulation 
depends strongly on the filling- the response of the double occupation is significantly different in 
the half-filled Mott insulator from the doped Fermi liquid region. 



PACS numbers: 03.75.Ss, 05.30.Fk, 34.50.-s, 71.10.Fd 

A number of key properties of strongly correlated elec- 
tron systems appear to be well described by simplified 
tight-binding Hamiltonians. For example, the square lat- 
tice Hubbard model, with one particle per site, is known 
to possess the long range antiferromagnetic order man- 
ifest in the parent compounds of high temperature su- 
perconductors, whose Cu02 sheets have square arrays of 
copper atoms with one hole per 3d shell. There are many 
analytic and numerical clues that suggest the doped Hub- 
bard model might also possess the d-wave superconduct- 
ing phase exhibited by the cuprates, as well as other non- 
trivial properties including stripes and pseudogap physics 
If this could be demonstrated rigorously, it would 
provide important insight into the mechanism of super- 
conductivity in these materials. 

Ultracold atomic systems offer an opportunity for 
closer connection between experiments and calculations 
for such model Hamiltonians. At present, experiments 
on fermionic atoms are exploring temperatures T which 
are of the order of the hopping integral Jq i probing cor- 
relations such as double occupancy, D, and short range 
spin order that develops at that temperature scale. In 
particular, the evolution of D with the ratio of interac- 
tion strength U to hopping Jg has been shown to in- 
dicate the presence of a Mott metal-insulator transition 
0, 01 • The presence of a Mott gap in the excitation spec- 
trum has also been inferred through peaks in D which 
arise through a dynamic modulation of the optical lat- 
tice depth V 

The possibility that such a modulation might provide 
a useful probe was first suggested by KoUath et al. [1], 
based on earlier work with bosonic systems Using a 
time dependent Density Matrix Renormalization Group 
method, it was shown that a peak existed in the induced 
double occupation at a frequency uj which matched the 
interaction strength U. In this treatment, the response 



kernel was approximated to include only changes 5 J in 
the hopping operator, neglecting corresponding variation 
SU in the on-site interactions. Within this approxima- 
tion, the authors emphasized that the measurement was 
sensitive to near neighbor spin correlations, and the ex- 
change gap, as well as the charge gap. 

This 'modulation spectroscopy' has been further ex- 
plored theoretically by HuberQ and SensarmaQ- In 
the former work, the frequency dependence of the shift 
in D was studied in the atomic and two particle limits, 
and within a slave boson mean field theory. The latter 
work focused on observing local antiferromagnetic order 
at the superexchange scale. As with the earlier study 
of Kollath, in both of these papers, the modulation was 
assumed to couple only to the kinetic energy. 

In this paper, we extend previous work by studying the 
effect of both the modulation of the tunneling strength 6 J 
and of the on-site interaction strength SU due to varying 
the optical lattice depth V, for the two dimensional re- 
pulsive fermionic Hubbard Hamiltonian. The modulation 
by SU is shown to be quite significant in the parameter 
range of interest to current experiments. We find that 
the filling of the system plays a very important role in 
the response. Crucially, through the use of Determinant 
Quantum Monte Carlo (DQMC) and the maximum 
entropy method [l^l , we provide results which treat 
the electron-electron correlations exactly. 

In the low energy limit, two species of repulsively inter- 
acting fermions confined to a periodic optical potential 
with wavelength A and amplitude V{t) can be described 
by the one-band Hubbard model 



H = -JK + UD- fiN, 



(1) 



where the hopping or kinetic-energy operator is K = 
E<ij>,<T[cLSff + h-c-], D = J2i %"-»4- is the double oc- 
cupancy, and N = fii-f + hii, the total number of par- 



2 



tides, with c\^{cia) the fermion creation (annihilation) 



operator, a =t,4. the spin index, hi^ = cl^Cj^, and the 
chemical potential. The hopping (J) and interaction (U) 
can be expressed as [lH J « (4 /v^) Ej^ v^l'^ exp(— 2-y/i;) 
and K, \yp2^(a,IX)ERV^I'^, where w = VjER is the 
ratio of lattice depth to recoil energy, and is the short 
ranged s-wave scattering length. 

It is clear from these expressions that a small time- 
dependent modulation of V changes both J and \J . Writ- 
ing V{€) ~ Vo + W smiijjt) and expanding J and U in 
the limit W < Vb yields H = Hq+ 6H sin(wt) with i?o 
given by Eq. ([T]) with J replaced by Jq and U hj Uq, and 
(5i? = —5JK + (5C/£) with the time-dependent perturba- 
tions 

5J = 

V4 V Vo 




For (5F > 0, we have 5J < and 5U > so that an 
increase in the optical lattice amplitude suppresses hop- 
ping and increases the Hubbard repulsion. We emphasize 
that one cannot a priori neglect 5 J or 5U as they can be 
of the same order of magnitude if using the experimental 
parameters as in Ref. [3- 

Our aim is to understand how such a simultaneous 
modulation of the hopping and interaction parameters, as 
provided by fcrmions in a time-dependent optical lattice, 
probes fermion correlations in the Hubbard model. To 
this end, we study the time dependence of the average 
double occupancy D{t) = (D). Within standard time- 
dependent perturbation theory, D{t) satisfies, to linear 
order. 



D{t) 



D{to)-i I dt' {[D {t),6H{t')]) sin ujt\ (3) 

'to 



iHot 



where {6)o = Z^^Tie-^^^O and d{t) = e^^^^Oc 
Equation ^ can be simplified by rewriting 5H in terms 
of as 5H = (5 J/ Jo) {Hq + Uo[a - l]D), with a = 

1 - ly^) • When inserted into Eq. the first 
term will give a vanishing contribution, leading to 

D{t) ^ D{ta) + -^{a^l) / dt' 6 J XDD{t - t') sin ut' , 

Jo Jto 

(4) 

where xoo(i-^') = -i{[d{t),d{t')])a 0{t -t'). Formally 
setting a = amounts to neglecting the modulation of 
the interaction term. In contrast, experimentally, a typ- 
ically varies within the range —0.41 < a < —0.28. The 
simplification leading to Eq. (|4]), can be generalized to 
show that XDD(i) = (Jo/C^o)'^XKK(i)i ^ f^'^t that we shall 
use below in our analysis. 

Numerically, we calculate the imaginary-time quantity 
Xdd{t) from Determinant Quantum Monte Carlo simu- 



lations [8| and analytically extrapolate to the correspond- 
ing imaginary part of the real frequency quantity Xdd('-^) 
by inverting 



du! 



(5) 



via the Maximum Entropy method [l^l- In Eq. [5] 
ivn = 2m:T is the bosonic Matsubara frequency, T is 
the temperature, and w the real frequency. 
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FIG. 1: (Color online). Top panel (a) shows data for half fill- 
ing, and panel (b) for a filling of n = 1.4, for a two-dimensional 
4x4 Hubbard lattice. Red curves (squares) show the quan- 
tity (1 — a)xuo, that appears in the linear response of the 
double occupancy, evaluated at zero Matsubara frequency as 
a function of Uo/Jo- Neglecting the modulation of the Hub- 
bard interaction amounts to setting a — 0, yielding a smaller 
result (black curve, circles). For comparison, the green dia- 
monds in the insert in both (a) and (b) are exact results for 
(1 — a)xDD for a two-site Hubbard model, a is determined 
by assuming Os/A — 0.0119, where Ua = 240ao (ao is Bohr 
radius) and A = 1, 064 nm (following Ref. Ql), thus a can be 
found as a single- valued function of Uo/Jo- 

To illustrate the importance of incorporating the mod- 
ulation of the interaction parameter J7, in Fig.[T]wc show 
the dependence with Uq / Jq of the double-occupancy re- 
sponse function XDD(ii'n = 0) (black curves), for n = 
{rii-f + nil) ^ ^-^ ^^'^ n ~ lA along with this quantity 
multiplied by (1 — a) (red curve). Therefore, the black 
curves is the result from modulating SJ only, while the 
red curve also includes the effect of modulating SU. The 
difference between the curves illustrates that SU should 
not be neglected. We observe from Fig. [T] that at half- 
filling (n=l), the double occupancy response is largest 
in the intermediate interaction region and decreases with 
increasing Uq/Jq. This is in striking contrast to the be- 
havior at n = 1.4, in which the double occupancy re- 
sponse is small at weak coupling and saturates at large 
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Uq/Jq. To confirm our numerical calculation, we analyt- 
ically solved the case of a two-site Hubbard model and 
found qualitatively similar behavior. (See green curves 
in Fig.m) 




FIG. 2: (Color online). The imaginary component of the 
double-occupancy susceptibility x'dd{'^)/N for Uo/Jo = 10.0, 
a 4 X 4 square lattice, and various values of inverse tempera- 
ture {P = 1/T). Panel (a) shows half-filling n = 1.0 results, 
and panel (b) a filling of n = 1.4. A'' = 16 is the system size. 

We now turn to the full frequency dependent dy- 
namical susceptibility, which determines the response to 
the dynamical modulation, showing its evolution as a 
function of temperature (expressed in terms of /3Jo = 
Jo/ksT) in Fig. [2j Panel (a) displays results at half- 
filling, where Mott-insulating physics dominates. At this 
filling the low frequency response is strongly suppressed 
for temperatures approaching zero (so that this quasi- 
peak represents thermally-excited states, not coherent 
excitations) , with the predominant response occurring at 
frequencies close to Uq. This energy scale, corresponding 
to the Mott gap, is consistent with recent experimen- 
tal results which find a strong response in the double 
occupancy when uj ^ Uq. The presence of the Mott gap 
also accounts for the much smaller values of x" in the top 
panels of Figs. [1] and [2] Panel (b) shows a filling n = 1.4, 
where an w = peak remains robust for T —> 0. We 
attribute this peak to the presence of gapless excitations 
reflecting Fermi liquid behavior in this region. The peak 
at high uj represents coherent excitations at the band-gap 
scale which should be the distance between the lower and 
upper Hubbard bands. 

In Fig. [3J we show the interaction dependence of 
Xdd(w). Panel (a) displays the half-filled case where the 
peaks are centered at Uq. In panel (b), filling n = 1.4, 
we include the case of a larger lattice size (6 x 6) to 
show that finite size effects are small. These results fur- 
ther verify the important role of filling in the response 
to dynamical modulation. Our findings can be quali- 
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FIG. 3; (Color online). Left column: The imaginary part of 
the double-occupancy susceptibility Xdd('^)/-^ for Uo/Jo = 
10 and 16. Panel (a) shows half-filling n = 1.0 results for a 
4x4 lattice, Uo/Jo = 10.0 (black solid curve) and Uo/Jo = 
16.0 (red solid curve). Panel (b) shows results for a filling 
n — 1.4 and for Uo/Jo = 10.0, 6x6 square lattice (orange 
curve), Uo/Jo = 10.0, 4x4 (green curve), and Uo/Jo ~ 16.0, 
4x4 (blue curve). Right column: The real-time double- 
occupancy response function ^dd (t) for a 4 x 4 square lattice 
at half filling (panel (c)) for Uo/Jo = 10.0 (black solid curve) 
and Uo/Jo ~ 16.0 (red solid); and for n = 1.4 (panel (d)) with 
Uo/Jo = 10.0 (green curve) and Uo/Jo = 16.0 (blue curve). 
All results are at a temperature T/Jo — 2/3. 



tatively reproduced by neglecting vertex corrections in 
Xkk and expressing the single particle Green's function 
in the Hubbard-I approximation. The latter corresponds 
to using a approximate self-energy of the form 



(6) 



Wc find that Xkk('^) (and hence XBui^)) possess poles 
at (jj ~ 0, ±A/(ek)'^ + 4C/q ?ia(l — "a), where ek is the en- 
ergy of a non-interacting quasiparticle with momentum 
k. In the low energy region, there are quasi-elastic peaks 
at approximately w ^ 0. Note that the peak vanishes at 
a; = because the imaginary part of the real frequency 
susceptibility is an odd function XkkI^'^) = ^Xkk('^)- 
In the high energy region, the peaks are located at 



roughly u! ^ Uo 



2Un 



Therefore, at half-filling, the 



peaks are at u) = Uq but they sit at higher frequencies 
away from half filling. 

We now turn to the question of how the features in 
Xdd(w) would be reflected in a experimental measure- 
ment of the double occupancy, by inserting our results 
for XDD(i) into Eq. (|4|). For this task, we need to ob- 
tain the real part of xdd(w) via Kramers-Kronig; upon 
Fourier transforming we find the real-time dynamical re- 
sponse functions for the double occupancy to be strik- 
ingly different at half filling and away from half filling, as 
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seen in panels (c) and (d) of Fig. [3l We see that filling 
n = 1 shows a response function that is tightly peaked at 
i -> 0, characterized by a single frequency scale oj ^ Uq, 
while at n = 1.4 we see a broad behavior dominated by 
the two distinct frequencies associated with u! ^ and 
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FIG. 4: (Color online). The frequency dependence of the 
double occupancy linear response for a 4 x 4 lattice, interac- 
tion strength Uo/Jo = 10.0, and temperature T/Jo = 2/3 . 
Panel (a) shows half-filling results; panel (b), n — 1.4. Solid 
(black) curves shows the amplitude D{lj) while the dashed 
(red) curves display the phase shift induced by the dy- 

namical modulation. 

As in standard linear response theory, the real and 
imaginary parts of xbi){^) correspond to the in-phase 
and out of phase parts of the response, respectively. 
Thus, to linear order, an oscillatory driving of the opti- 
cal lattice potential yields an oscillatory response at the 
same frequency, but with a phase lag characterized by 
the ratio of tan0(a;) = XuDi^) / XDoiiiil- This response 



has recently been observed directly [12[. We can then 



(7) 



write the time-dependent double occupancy as 

D{t) = D{0) + D{uj) sin[ujt - 0(a;)], 

where D{uj) = Uo/Jo{a - 1)W|xdd(w)|. We plot D{uj) 
and in Fig. |4]for the case of Uq/Jq = 10. We first 
note that, at low frequency w 0, Eq. ([7]) implies the 
time dependence of D{t) to be precisely tt out of phase 
with 6V{t). Therefore, an adiabatic increase of the opti- 
cal lattice amplitude leads to a corresponding suppression 
of the double occupancy. At higher uj these plots show 
how the time-dependent linear response of the double oc- 
cupancy probes the underlying fermion correlations. As 
we expected, the half filled case shows the strongest re- 
sponse when the driving frequency lu ^ U, and with a 
phase that is shifted, by « 7r/2, relative to the imposed 
modulation. At (n) = 1.4, however, the predominant re- 
sponse is for = 0, with phase shift « 0. 



In conclusion, we have investigated the dynamical 
properties of fermions in an optical lattice, realized by 
the Hubbard model subject to a periodic optical lattice 
modulation. We show that the modulation of the on- 
site interaction cannot be neglected and that, even at 
the level of linear response, the dynamical double occu- 
pancy provides a sensitive probe of fermion correlations. 
Recent cold-atom experiments studying the dynami- 
cal modulation of the optical lattice find a linear in time 
contribution to the double occupancy, known to emerge 
at quadratic order in the modulation parameter SV [J]. 
Thus, we expect that our linear-response results apply 
at smaller SV/Vq, or after subtracting off this ^-linear 
contribution to focus on the oscillatory component. A 
future extension of our work will analyze the linear and 
quadratic-order contributions in detail. In addition, the 
effects of inhomogeneity due to trapping effects is an issue 
for future calculations. 
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